    
    DataGen = function(n, SEED=0){
        
        # set.seed(SEED)
        # set.seed(seeds[seed.index]) #DEBUG
        
        l0.1        = rep(1,n)  
        l0.2        = rep(1,n) 
        l0.2.tilde  = rep(1,n) 
        if (l0.type == "continuous"){
          l0.2        = runif(n,-2,2)
        }else if (l0.type == "discrete"){
          l0.2        = rbinom(n,1,0.5)
          l0.2.tilde       = rbinom(n,1,0.5)
          # if (gop.type == "incorrect"){
          #     l0.2.tilde       = rbinom(n,1,0.5)
          # }
        }
        l0          = cbind(l0.1,l0.2) 
        l0.tilde    = cbind(l0.1,l0.2.tilde) 

        p.a0.true   = expit(l0 %*% gamma.true["a0",1:2])
        a0          = rbinom(n,1,p.a0.true)
        p.l1.true   = numeric(n)
        p.l1.true[a0==1] = expit(l0[a0==1,] %*% beta.true[4,])
        p.l1.true[a0==0] = expit(l0[a0==0,] %*% beta.true[5,])
        l1          = rbinom(n,1,p.l1.true)
        p.a1.true   = expit(cbind(l1,a0,l0) %*% gamma.true["a1",])
        a1          = rbinom(n,1,p.a1.true)
        
        theta.true  = exp(l0 %*% t(alpha.true))
        phi123.true = exp(l0 %*% t(beta.true[1:3,]))
        phi45.true  = expit(l0 %*% t(beta.true[4:5,]))
        phi.true    = cbind(phi123.true, phi45.true)
        if (vectorize.get.prob){
          p.y.cond.true = GetProb(theta.true,
                                  phi.true)
        }else{
          p.y.cond.true = t(sapply(1:n, function(i) GetProb(theta.true[i,],
                                                            phi.true[i,])))
        }
        colnames(p.y.cond.true) = c("p000","p001","p010","p011","p100",
                                    "p101","p110","p111")
        a0l1a1 = 4*a0 + 2*l1 + a1 + 1
        # table(a0l1a1) #DEBUG
        p.y.true = sapply(1:n, function(i) p.y.cond.true[i,a0l1a1[i]])
        y      = rbinom(n,1,p.y.true)
        
        if (l0.type == "continuous"){
          data <- cbind(l0, a0,l1,a1,y)
          cnt <- rep(1,n)
        }else if (l0.type == "discrete"){
          data.raw = list(l0.1 = 1, l0.2 = c(0,1), a0 = c(0,1), l1 = c(0,1), 
                          a1 = c(0,1), y = c(0,1))
          data           = ParaValues(data.raw)
          dim(data)      = c(32, 6)
          colnames(data) = names(data.raw)
          cnt            = as.numeric(table(l0.1, l0.2, a0, l1, a1, y))
          
          if (gop.type == "incorrect"){
              data.raw = list(l0.1 = 1, l0.2 = c(0,1),l0.2.tilde = c(0,1),
                              a0 = c(0,1), l1 = c(0,1), 
                              a1 = c(0,1), y = c(0,1))
              data           = ParaValues(data.raw)
              dim(data)      = c(64, 7)
              colnames(data) = names(data.raw)
              cnt            = as.numeric(table(l0.1, l0.2,l0.2.tilde, a0, l1, a1, y))
          }
        }
        
        return(list(data = data, cnt = cnt, 
                    a1.vec = a1, a1.cov = cbind(l1,a0,l0),
                    a0.vec = a0, a0.cov = l0, 
                    l0.2.vec=l0.2, l1.vec = l1,
                    y.vec = y, y.cov=cbind(l1,a0,l0,a1),
                    y.cov.tilde=cbind(l1,a0,l0.tilde,a1),
                    p.a1.true=p.a1.true,
                    p.y.cond.true=p.y.cond.true,p.a0.true=p.a0.true,
                    p.l1.true=p.l1.true))
    }
    